home *** CD-ROM | disk | FTP | other *** search
/ MacFormat 1995 April / macformat-023.iso / Shareware City / Developers / RLaB / toolbox / rk4.r < prev    next >
Encoding:
Text File  |  1994-02-21  |  1.4 KB  |  78 lines  |  [TEXT/RLAB]
  1. //
  2. //  RLaB 4th Order Runge-Kutta Integrator
  3. //  Fixed step size.
  4. //  Ian Searle, 6/22/92
  5.  
  6. //
  7. //  Inputs to rk4()
  8. //  eval:    User-Function to evaluate state-derivative.
  9. //        has arguments (time, y)
  10. //  tzero:        Start time
  11. //  tfinal:       Finish time
  12. //  y0:              Initial conditions
  13. //  dt:              Integration time step
  14.  
  15. //  Output:       Matrix of time, and state vector.
  16. //
  17.  
  18. static (collapse, rrk4);
  19.  
  20. rk4 = function(eval, tzero, tfinal, y0, dt) 
  21. {
  22.   local(i, y, lout, I);
  23.  
  24.   y = y0[:];
  25.  
  26.   lout = << [tzero, y0.'] >>; I = 2;
  27.   for(i in tzero:(tfinal-dt):dt) 
  28.   {
  29.     lout.[I] = [i+dt, (y = rrk4 (y, i, dt, eval))'];
  30.     I++;
  31.   }
  32.   return collapse (lout);
  33. };
  34.  
  35. rrk4 = function(y, x, h, eval) 
  36. {
  37.   local(i, xh, hh, h6, dydx, dym, dyt, yt, yout);
  38.  
  39.   // Initialize
  40.   hh = h*0.5;
  41.   h6 = h/6;
  42.   xh = x + hh;
  43.  
  44.   dydx = eval(x, y);      // 1st step
  45.   yt = y + hh*dydx;
  46.  
  47.   dyt = eval(xh, yt);     // 2nd step
  48.   yt = y + hh*dyt;
  49.  
  50.   dym = eval(xh, yt);     // 3rd step
  51.   yt = y + h*dym;
  52.   dym = dym + dyt;
  53.  
  54.   dyt = eval(x+h, yt);    // 4th step
  55.   yout = y + h6*(dydx + dyt + 2.0*dym);
  56.  
  57.   return yout;
  58. };
  59.  
  60. //
  61. // Collapse a LIST of matrices into a single matrix.
  62. //
  63.  
  64. collapse = function ( LIST )
  65. {
  66.   local (i, m, row, col);
  67.  
  68.   row = size (LIST.[members (LIST)[1]])[1];
  69.   col = size (LIST.[members (LIST)[1]])[2];
  70.   m = zeros (length (LIST)*row, col);
  71.  
  72.   for (i in 1:length (LIST))
  73.   {
  74.     m[i;] = LIST.[i];
  75.   }
  76.   return m;
  77. };
  78.